# Third tutorial¶

## Crystalline silicon.¶

This tutorial aims at showing you how to get the following physical properties, for an insulator:

• the total energy
• the lattice parameter
• the band structure (actually, the Kohn-Sham band structure)

You will learn about the use of k-points, as well as the smearing of the plane-wave kinetic energy cut-off.

Important

All the necessary input files to run the examples can be found in the ~abinit/tests/ directory where ~abinit is the absolute path of the abinit top-level directory.

To execute the tutorials, you are supposed to create a working directory (Work*) and copy there the input files and the files file of the lesson.

The files file ending with _x (e.g. tbase1_x.files) must be edited every time you start to use a new input file. You will discover more about the files file in section 1.1 of the help file.

To make things easier, we suggest to define some handy environment variables by executing the following lines in the terminal:

export ABI_HOME=Replace_with_the_absolute_path_to_the_abinit_top_level_dir

export ABI_TESTS=$ABI_HOME/tests/ export ABI_TUTORIAL=$ABI_TESTS/tutorial/           # Files for base1-2-3-4, GW ...
export ABI_TUTORESPFN=$ABI_TESTS/tutorespfn/ # Files specific to DFPT tutorials. export ABI_TUTOPARAL=$ABI_TESTS/tutoparal/         # Tutorials about parallel version
export ABI_TUTOPLUGS=$ABI_TESTS/tutoplugs/ # Examples using external libraries. export ABI_PSPDIR=$ABI_TESTS/Psps_for_tests/       # Pseudos used in examples.

export PATH=$ABI_HOME/src/98_main/:$PATH


The examples in this tutorial will use these shell variables so that one can easily copy and paste the code snippets into the terminal (remember to set ABI_HOME first!)

The last line adds the directory containing the executables to your PATH so that one can invoke the code by simply typing abinit in the terminal instead of providing the absolute path.

Finally, to run the examples in parallel with e.g. 2 MPI processes, use mpirun (mpiexec) and the syntax:

mpirun -n 2 abinit < files_file > log 2> err


The standard output of the application is redirected to log while err collects the standard error (runtime error messages, if any, are written here).

This tutorial should take about 1 hour.

Visualisation tools are NOT covered in this tutorial. Powerful visualisation procedures have been developed in the Abipy context, relying on matplotlib. See the README of Abipy and the Abipy tutorials.

## Computing the total energy of silicon at a fixed number of k-points¶

Before beginning, you might consider working in a different subdirectory, as for tutorial 1 or 2. Why not Work3?

The file tbase3_x.files lists the file names and root names. You can copy it in the Work3 directory and change it as you did in the first and second tutorials. You can also copy the file tbase3_1.in inside the Work3 directory with:

while $ABI_TUTORIAL/Refs/tbase3_5.out is a reference output file. You should find the band structure starting at (second dataset):  Eigenvalues ( eV ) for nkpt= 40 k points: kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord) -3.78815 -1.15872 4.69668 4.69668 7.38795 9.23867 9.23867 13.45707 kpt# 2, nband= 8, wtk= 1.00000, kpt= 0.4500 0.0000 0.0000 (reduced coord) -3.92759 -0.95774 4.71292 4.71292 7.40692 9.25561 9.25561 13.48927 kpt# 3, nband= 8, wtk= 1.00000, kpt= 0.4000 0.0000 0.0000 (reduced coord) -4.25432 -0.44393 4.76726 4.76726 7.46846 9.31193 9.31193 13.57737 kpt# 4, nband= 8, wtk= 1.00000, kpt= 0.3500 0.0000 0.0000 (reduced coord) -4.64019 0.24941 4.85732 4.85732 7.56855 9.38323 9.38323 13.64601 ....  One needs a graphical tool to represent all these data. In a separate file (_EIG), you will find the list of k-points and the eigenenergies (the input variable prteig is set by default to 1). Even without a graphical tool we will have a quick look at the values at L, $\Gamma$, X and $\Gamma$ again:  kpt# 1, nband= 8, wtk= 1.00000, kpt= 0.5000 0.0000 0.0000 (reduced coord) -3.78815 -1.15872 4.69668 4.69668 7.38795 9.23867 9.23867 13.45707 kpt# 11, nband= 8, wtk= 1.00000, kpt= 0.0000 0.0000 0.0000 (reduced coord) -6.17005 5.91814 5.91814 5.91814 8.44836 8.44836 8.44836 9.17755 kpt# 23, nband= 8, wtk= 1.00000, kpt= 0.0000 0.5000 0.5000 (reduced coord) -1.96393 -1.96393 3.00569 3.00569 6.51173 6.51173 15.95524 15.95524 kpt# 40, nband= 8, wtk= 1.00000, kpt= 1.0000 1.0000 1.0000 (reduced coord) -6.17005 5.91814 5.91814 5.91814 8.44836 8.44836 8.44836 9.17755  The last $\Gamma$ is exactly equivalent to the first $\Gamma$. It can be checked that the top of the valence band is obtained at $\Gamma$ (=5.91814 eV). The width of the valence band is 12.09 eV, the lowest unoccupied state at X is 0.594 eV higher than the top of the valence band, at $\Gamma$. The Si is described as an indirect band gap material (this is correct), with a band-gap of about 0.594 eV (this is quantitatively quite wrong: the experimental value 1.17 eV is at 25 degree Celsius). The minimum of the conduction band is even slightly displaced with respect to X, see kpt # 21. This underestimation of the band gap is well-known (the famous DFT band-gap problem). In order to obtain correct band gaps, you need to go beyond the Kohn-Sham Density Functional Theory: use the GW approximation. This is described in the first tutorial on GW. For experimental data and band structure representation, see the book by M.L. Cohen and J.R. Chelikowski [Cohen1988]. Important There is a subtlety that is worth to comment about. In non-self-consistent calculations, like those performed in the present band structure calculation, with iscf = -2, not all bands are converged within the tolerance tolwfr. Indeed, the two upper bands (by default) have not been taken into account to apply this convergence criterion: they constitute a buffer. The number of such buffer bands is governed by the input variable nbdbuf. It can happen that the highest (or two highest) band(s), if not separated by a gap from non-treated bands, can exhibit a very slow convergence rate. This buffer allows achieving convergence of important, non-buffer bands. In the present case, 6 bands have been converged with a residual better than tolwfr, while the two upper bands are less converged (still sufficiently for graphical representation of the band structure). In order to achieve the same convergence for all 8 bands, it is advised to use nband=10 (that is, 8 + 2). ## Using AbiPy to automate the most boring steps¶ The AbiPy package provides several tools to facilitate the preparation of band structure calculations and the analysis of the output results. First of all, one can use the abistruct script with the kpath command to determine a high-symmetry k-path from any file containing structural information (abinit input file, netcdf output files etc.). The high-symmetry k-path follows the conventions described in [Setyawan2010]. Let’s try with: abistruct.py kpath tbase3_5.in # Abinit Structure natom 2 ntypat 1 typat 1 1 znucl 14 xred 0.0000000000 0.0000000000 0.0000000000 0.2500000000 0.2500000000 0.2500000000 acell 1.0 1.0 1.0 rprim 0.0000000000 5.1080000000 5.1080000000 5.1080000000 0.0000000000 5.1080000000 5.1080000000 5.1080000000 0.0000000000 # K-path in reduced coordinates: # tolwfr 1e-20 iscf -2 getden ?? ndivsm 10 kptopt -11 kptbounds +0.00000 +0.00000 +0.00000 #$\Gamma$+0.50000 +0.00000 +0.50000 # X +0.50000 +0.25000 +0.75000 # W +0.37500 +0.37500 +0.75000 # K +0.00000 +0.00000 +0.00000 #$\Gamma\$
+0.50000  +0.50000  +0.50000 # L
+0.62500  +0.25000  +0.62500 # U
+0.50000  +0.25000  +0.75000 # W
+0.50000  +0.50000  +0.50000 # L
+0.37500  +0.37500  +0.75000 # K
+0.62500  +0.25000  +0.62500 # U
+0.50000  +0.00000  +0.50000 # X


To visualize the band structure stored in the GSR.nc file, use the abiopen script and the command line:

abiopen.py tbase3_5o_DS2_GSR.nc --expose -sns=talk


It is also possible to compare multiple GSR files with the abicomp script and the syntax

abicomp.py gsr tbase3_5o_DS1_GSR.nc tbase3_5o_DS2_GSR.nc -e -sns=talk


to produce the following figures:

For further details about the AbiPy API and the GSR file, please consult the GsrFile notebook .